clear all;
%% import data file
%rawdata=importdata('D:\yaomingxi\Desktop\2012-03-10\ac585biosmcc_100mMtris8_150mMnacl_pcatrolox_bead04_recalibrated_fluctuations.txt');

% For large data
% This code loads in the contents of the data file with textscan a block at a
% time and converts the cellarray from textscan into a double array

filename='G:\unfolding_vinculin\ntasur_m270strept_bead120_200fps_s1_1xpbs_constforce.txt'; 
FormatString='%f %f %f %f %f';

% Parameters
Headers=2;
NumLines=2607526; % Total number of lines to read
BlockLines=10000; % Size of block
NumBlocks=NumLines/BlockLines; % Number of blocks

% Open file
fid = fopen(filename); % Open file

% Preallocate data
mydata=zeros(NumLines,5); % Pre-allocate space for data

% Get Headers
cellchunk=textscan(fid,'%s',2,'delimiter','\n'); % Get first 3 lines

% Read in blocks
for Block=1:NumBlocks
    mydata((Block-1)*BlockLines+(1:BlockLines),:) = ... 
	cell2mat(textscan(fid,FormatString,BlockLines,'delimiter','\t'));
    disp(['Block ' num2str(Block) ]); % Display current block to show progress
end

% Close file
fclose(fid);


%get individual columns
time=mydata(:,1);
magnet=mydata(:,2);
extension=mydata(:,3);
dx=mydata(:,4);
dy=mydata(:,5);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Parameters
position_overstretching=12.45;
%c=65/(exp(-(13.5-position_overstretching)/0.36)+0.48*exp(-(13.5-position_overstretching)/1.12));
c=280;
span = 20; % Size of the averaging smooth window
%size of the averaging window
window_size=20;
%threshold step size in nm
threshold=8;


%% plot timelapse curve
figure
subplot(2,4,1)
plot(time,extension,'b-');
hold on

window = ones(span,1)/span; 
smoothed_extension = convn(extension,window,'same');
plot(time,smoothed_extension,'r-');
legend('Data','Smoothed Data')
title('Time Lapse Extension');
xlabel('Time(s)');
ylabel('Extension(nm)');
hold off

%find when magnet were changed get segments of individual forces for
%processing
data_length=length(time);
flag_changeMagPosition=1;
position_magnetStart=-999;
MagnetPositions=zeros(0,4);
for i=1:data_length
    position_magnetCurrent=magnet(i);
    if position_magnetCurrent==position_magnetStart
        flag_changeMagPosition=0;
        count_magPosition=count_magPosition+1;
    else
        flag_changeMagPosition=1;
    end
    
    if flag_changeMagPosition==1
        if i~=1
            force_current=c*(exp(-(13.5-position_magnetStart)/0.36)+0.48*exp(-(13.5-position_magnetStart)/1.12));
            MagnetPositions=[MagnetPositions; [position_magnetStart point_positionStart count_magPosition force_current]];
        end
        point_positionStart=i;
        position_magnetStart=magnet(i);
        count_magPosition=0;    
    end
    if i==data_length
            force_current=c*(exp(-(13.5-position_magnetStart)/0.36)+0.48*exp(-(13.5-position_magnetStart)/1.12));
            MagnetPositions=[MagnetPositions; [position_magnetStart point_positionStart count_magPosition force_current]];
    end
end

%Get statistics at individual forces 
MatrixStatistics=zeros(0,8);
for i=1:length(MagnetPositions)
    extension_current=extension(MagnetPositions(i,2):MagnetPositions(i,2)+MagnetPositions(i,3));
    dx_current=dx(MagnetPositions(i,2):MagnetPositions(i,2)+MagnetPositions(i,3));
    dy_current=dy(MagnetPositions(i,2):MagnetPositions(i,2)+MagnetPositions(i,3));
    mean_extension_current=mean(extension_current);
    std_extension_current=std(extension_current);
    mean_dx_current=mean(dx_current);
    std_dx_current=std(dx_current);
    mean_dy_current=mean(dy_current);
    std_dy_current=std(dy_current);
    MatrixStatistics=[MatrixStatistics; [MagnetPositions(i,1) MagnetPositions(i,4) mean_extension_current std_extension_current mean_dx_current std_dx_current mean_dy_current std_dy_current]];
end

%Plot the force extension curve and force variance curve

%identify differet passes and assign them different color

colorprofile=zeros(length(MatrixStatistics),3);
color_current=1;
pass_current=1;
flag_peak=0;
position_peaks=1;
minpeakamp=9.5;

for i=1:length(MatrixStatistics)
    if (MatrixStatistics(i,1)>minpeakamp) && (flag_peak==0)
        position_peak_start=i;
        flag_peak=1;
    else if (MatrixStatistics(i,1)<minpeakamp) && (flag_peak==1)
            position_peak_end=i;
            flag_peak=0;
            peak_value=max(MatrixStatistics(position_peak_start:position_peak_end,1));
            for j=position_peak_start:position_peak_end
                if MatrixStatistics(j,1)==peak_value
                    tmp_peaks=j;
                end
            end
            if position_peaks(end)~=1
                low_peak_value=min(MatrixStatistics(position_peaks(end):tmp_peaks));
                for j=position_peaks(end):tmp_peaks
                    if MatrixStatistics(j,1)==low_peak_value
                        low_tmp_peaks=j;
                    end
                end
                position_peaks=[position_peaks;low_tmp_peaks];

            end
            position_peaks=[position_peaks;tmp_peaks];
            
        end
    end
end
colorlist=jet(length(position_peaks));
for i=1:(length(position_peaks)-1)
    colorprofile(position_peaks(i):position_peaks(i+1),:)=repmat(colorlist(color_current,:),length(colorprofile(position_peaks(i):position_peaks(i+1),:)),1);
    color_current=color_current+1;
end
colorprofile(position_peaks(end):end,:)=repmat(colorlist(color_current,:),length(colorprofile(position_peaks(end):end,:)),1);


%plot force fluctuation and force extension curves
subplot(2,4,2)
scatter(MatrixStatistics(:,2),MatrixStatistics(:,3),10,colorprofile,'filled')
title('Force vs Extension');
xlabel('Forces (pN)');
ylabel('Extension (nm)');
subplot(2,4,3)
scatter(MatrixStatistics(:,2),MatrixStatistics(:,4),10,colorprofile,'filled')
title('Force vs stdZ');
xlabel('Forces (pN)');
ylabel('stdZ (nm)');
subplot(2,4,4)
scatter(MatrixStatistics(:,2),MatrixStatistics(:,5),10,colorprofile,'filled')
title('Force vs dx');
xlabel('Forces (pN)');
ylabel('mean dx (nm)');
subplot(2,4,5)
scatter(MatrixStatistics(:,2),MatrixStatistics(:,6),10,colorprofile,'filled')
title('Force vs stdX');
xlabel('Forces (pN)');
ylabel('stdX (nm)');
subplot(2,4,6)
scatter(MatrixStatistics(:,2),MatrixStatistics(:,7),10,colorprofile,'filled')
title('Force vs dy');
xlabel('Forces (pN)');
ylabel('mean dy (nm)');
subplot(2,4,7)
scatter(MatrixStatistics(:,2),MatrixStatistics(:,8),10,colorprofile,'filled')
title('Force vs stdY');
xlabel('Forces (pN)');
ylabel('stdY (nm)');

% subplot(2,4,8)
% scatter(MatrixStatistics(5,3),MatrixStatistics(5,5),10,colorprofile,'filled')
% title('Extension vs dx');
% xlabel('Forces (pN)');
% ylabel('stdY (nm)');
%checking the unfolding signals find the points where the mean unfolding
%peaks are over certain point



%% Analyze steps
flag_step=0;
step_list=zeros(0,4);
unfold_list=zeros(0,4);
refold_list=zeros(0,4);
step_diff_max=0;
for i=1:length(MagnetPositions)
    extension_current=extension(MagnetPositions(i,2):MagnetPositions(i,2)+MagnetPositions(i,3));
    for j=(window_size+1):(length(extension_current)-window_size)
        previous=mean(extension_current((j-window_size):j));
        later=mean(extension_current(j:(j+window_size)));
        step_diff_current=later-previous;
        if abs(step_diff_current)>threshold 
                flag_step=1;
        end
        if flag_step==1
            if abs(step_diff_current)<=threshold 
                flag_step=0;
                step_list=[step_list;[MagnetPositions(i,1) MagnetPositions(i,4) step_position step_diff_max]];
                if step_diff_max>0
                    unfold_list=[unfold_list;[MagnetPositions(i,1) MagnetPositions(i,4) step_position step_diff_max]];
                else
                    refold_list=[refold_list;[MagnetPositions(i,1) MagnetPositions(i,4) step_position abs(step_diff_max)]];
                end
                step_diff_max=0;
            else
                if abs(step_diff_current)>abs(step_diff_max)
                    step_diff_max=step_diff_current;
                    step_position=MagnetPositions(i,2)+j-1;
                end
            end
            
        end
    end
end

%% plot histogram of the stepsizes

%define bin sizes by indicating the center for each bin
histcenters={-30:2:30 0:0.2:20};

X=[step_list(:,4) step_list(:,2)];
UNfold=[unfold_list(:,4) unfold_list(:,2)];
REfold=[refold_list(:,4) refold_list(:,2)];
figure
subplot(2,2,1)
hist3(X,histcenters)
title('Steps');
xlabel('Step size'); ylabel('Force');
set(gcf,'renderer','opengl');
set(get(gca,'child'),'FaceColor','interp','CDataMode','auto');

subplot(2,2,2)
hist3(UNfold,histcenters)
title('Unfold Steps');
xlabel('Step size'); ylabel('Force');
set(gcf,'renderer','opengl');
set(get(gca,'child'),'FaceColor','interp','CDataMode','auto');

subplot(2,2,3)
hist3(REfold,histcenters)
title('REfold Steps');
xlabel('Step size'); ylabel('Force');
set(gcf,'renderer','opengl');
set(get(gca,'child'),'FaceColor','interp','CDataMode','auto');

%get extension histograms of the forces that have steps
%force_with_steps=unique(step_list(:,1));

%% plot hisgram for the forces with high fluctuation

xbins=170:0.5:260;



for i=1:3:length(MagnetPositions)
    hFig=figure();
    set(hFig,'Visible','off');
    subplot(3,2,1)
    plot(time(MagnetPositions(i,2):(MagnetPositions(i,2)+MagnetPositions(i,3))),extension(MagnetPositions(i,2):(MagnetPositions(i,2)+MagnetPositions(i,3))));
    subplot(3,2,2)
    histnorm(extension(MagnetPositions(i,2):(MagnetPositions(i,2)+MagnetPositions(i,3))),xbins)
    if (i+1)<=length(MagnetPositions)
        subplot(3,2,3)
        plot(time(MagnetPositions(i+1,2):(MagnetPositions(i+1,2)+MagnetPositions(i+1,3))),extension(MagnetPositions(i+1,2):(MagnetPositions(i+1,2)+MagnetPositions(i+1,3))));
        subplot(3,2,4)
        histnorm(extension(MagnetPositions(i+1,2):(MagnetPositions(i+1,2)+MagnetPositions(i+1,3))),xbins)
    end
    if (i+2)<=length(MagnetPositions)
        subplot(3,2,5)
        plot(time(MagnetPositions(i+2,2):(MagnetPositions(i+2,2)+MagnetPositions(i+2,3))),extension(MagnetPositions(i+2,2):(MagnetPositions(i+2,2)+MagnetPositions(i+2,3))));
        subplot(3,2,6)
        histnorm(extension(MagnetPositions(i+2,2):(MagnetPositions(i+2,2)+MagnetPositions(i+2,3))),xbins)
    end
    print(hFig,'-dpsc2','-append','histogram.ps');
end



    
    
